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Abstract 

We consider a resonantly perturbed system of coupled nonlinear 
oscillators with small dissipation and outer periodic perturbation. We 
show that for the large time t ~ one component of the system is de- 
scribed in the main by the inhomogeneous Mathieu equation while the 
other component represents pulsation of large amplitude. A Hamilto- 
nian system is obtained which describes the behaviour of the envelope 
in the main. The analytic results agree to numerical simulations. 



^ '■ 1 Preliminaries 

in 

■ The paper deals with a system of coupled oscillators. This system has a 

large number of applications. In the framework of classical theory it can 
be used e.g. for describing oscillations of electrons and atoms under forced 
O ! combinative scattering, see 

For instance consider isotropic medium and denote by u the normal co- 
ordinate of oscillations of electrons in atoms and by v the normal coordinate 
^ ■ of oscillations of atoms in molecules. The frequency of oscillations u cor- 

responds to the ultraviolet scale range while the frequency of oscillations v 
lies in the ultrared scale range. To a linear approximation the oscillations of 
electrons and atoms occur independently of each other. 

In the general case on taking third order terms into account the expression 
for the potential energy has the form 

U{u, v) = ^au^ + ^f3v^ + /3iv^ + /32uv^ + /^sw^ + ^Wv, 



*Ufa State Petroleum Technical University (glebskie@gmail.com) 
■^Institute of Mathematics USC RAS (ok@ufanet.ru) 

■^Institute of Mathematics, Potsdam University (tarkhanov@math.uni-potsdam.de) 



1 



where a and /3 are constraint elasticity coefficients in a molecule, the coef- 
ficient j3i determines the nonlinearity of quasielasticity of oscillations f , (32 
determines the parametric stimulation of oscillations v by means of electron 
oscillations, /Ss determines the nonlinearity of electron oscillations which are 
responsible for the generation of optical harmonics. The coefficient 7 corre- 
sponds to the nonlinear interaction of electrons and nucleus and it determines 
the process of combinative scattering. Here we will consider the process of 
forced combinative scattering which is related to the presence of antihar- 
monic terms in the expression for the potential energy U . We assume in the 
sequel that the constants a, /3, /32, are essentially less than 7. 

Under the assumptions that both u and v are small we get the system 

u" + UiU + uj'^u 
v" + U2V' + f2^v 

see [3]. Here m and M are reduced masses of electronic and atomic oscilla- 
tors, 00 = \foiJrn and i? = \f^jM are eigen oscillation frequences, 7 is the 
constant that describes the nonlinear interaction of electrons and nucleus, 
e is the electron charge, E stands for the electric field of light waves, and 
1^1, V2 are dissipation coefficients. On assuming the atomic dissipation to be 
essentially less that the electronic dissipation, we will actually neglect V2 in 
the sequel. 

Change the dependent variables in system fll.ip by x = —u and y = —v. 
On taking the radiation damping in an electronic oscillator into account, see 
[3], [5], we arrive at the system 

+ ux'^ + u'^x 

y'u + ^'y 

where e = 7/m is a small positive parameter, 6 = m/(2M), the quantity 
{e/m)E = — A cos (i?t/2) is periodic in t, and u, i?. A, v are positive con- 
stants. 

Other examples of systems of coupled oscillators concern the description 
of propagation of surface and interior gravitational waves, see [6], [7]. 

The simplest analogue of system fll.ip is a nonlinear oscillator with res- 
onance pumping. The behaviour of such an oscillator is characterised at 
large times by beginning of permanent periodic nonlinear oscillations of the 
envelope whose amplitude is large when compared with the pumping. This 
process is described by the primary resonance equation, see for instance [8], 
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Fig. 1: Linear growth of the ?/ -component of solution of (11.11) in the initial 
stage. The solution amplitude is well approximated by the straight line 
y = ^eAHl{^u? - n'^f. 

The system fll.ip belongs to the other class of resonantly perturbed prob- 
lems. The pumping fails here to occur directly, i.e. by including the pertur- 
bation of eigen frequency to the right-hand side of the equation, but rather 
by means of nonlinear interaction of oscillators. In this situation the Mathieu 
equation proves to be of crucial role. 

In the paper an analogue of the nonlinear resonance equation for the 
envelope of resonant component is derived in the case of special nonlinear 
coupling in system fll.ip . 

In the initial stage the solution amplitude growths linearly, see Fig. [H 
This growth corresponds to the linear resonance in the second equation of 
the system. The envelope is well approximated by the straight line y = 
'ieA^t/ [Au"^ — Q'^Y. For large values of y the nonlinear effects of interaction 
between oscillators become essential. 

In Fig. [2] the results of numerical simulation of the behaviour of system 
(ll.ip of weakly coupled oscillators are shown under zero initial data. The 
parameter values are e = 0.2, i7 = l, ci; = 3, ^4 = 1 and v = 0.1 Numerical 
simulations are realized by the Adams method. 

The main result of the paper consists in describing dynamics of the en- 
velope of oscillations for large times. It turns actually out that the envelope 
oscillates as well. We show formulas for the amplitude and period of these 
oscillations. 

The paper consists of 5 sections. In Section [2] we formulate the results 
obtained. In Section [3] an asymptotic solution to the initial system fll.ip is 
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Fig. 2: Numerical simulation of the behaviour of system (11. ip with e = 
0.2, i? = l,uj = 3, A = 1 and u = 0.1. One can see pulsations of the y- 
component of solution of (11.11) . In the initial stage the solution amplitude 
is well approximated by the straight line y = AeA^t/ {Au'^ — Q'^Y. For large 
values of y the nonlinear effects of interaction between oscillators become 
essential. 

constructed. In Section H] we study in detail the physically interesting case 
where the frequency quotient cj/i? ^ 1 is large enough. In this case explicit 
formulas for the amplitude and period of oscillations are given. In the last 
Section [5] we shortly comment on the results. 

2 Results 

The asymptotic solution of system (11.11) for the large time t = 0(£:~^) has 
the form 

x{t,e) ~ xo{t,e'^t), 
y{t,e) ~ e-'Y{t,eH) 

for 5 — )■ 0, where e = — and 

Y{t, eh) = 2 3? k{eh) exp(zf?t), 

standing for the real part of a complex number z. The function xq satisfies 
the inhomogeneous Mathieu equation 

dl^XQ + A^dsXo + 4(A^ - \k\ cos2s)a;o = 4-^ cos (s - ^ arg/cV (2.1) 
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where 

f]t + arg k u 

' = 2 ' ^ = 7?- 

Write k = ki + tk2, where ki = 3ft/c and k2 = '^k are the real and 
imaginary parts of the function k, respectively. They are solutions of the 
averaged system 

drki = — lim — - / x1{'d,T) smd d-d, 

T-foo T Jo 

drk2 = lim — - / Xq('^, r) COS'!? (ii?, 

T^oo T Jo 

where r = e'^t. 

In the case A = w/i? ^ 1 which is important for applications the last 
system simplifies and reduces to the Hamiltonian system with Hamiltonian 
of form 



and i^i = ^,K2 = ^. 

The maximal value of the amplitude y for H = Hq and for t = O^s""^) 
with e -> is evaluated by 

max \y\ ~ — max{|r_|, \r^\}, 

where A ^ 1 and 



~ AtiQ^Hq ■ 

The period is in turn determined by the formula 



the integration being over the cycle c in the complex plane of the variable 
K = K^ + 1K2 given by H{K^, K2) = Hq. 
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3 Formal constructions for e <^ I 

In this section we write down the problems for determining the coefficients of 
asymptotic solutions to system (11 .ip . For convenience, we introduce a new 
independent variable 9 = Qt. On designating 

\ - ^ - ^ f - ^ - ^ 

7?' ■'~7p' ^~72' 

we get 

Xgg + fIXg + A^X = S Xy + f COS f , , 

y'ee + y = edx"^. 

The solution of this system is constructed by the two scales method (see 
[9]) in the form 

x{e,T,e) = Xo{9,t), 

y{9,T,e) = e-'Y{9,r)+ey^{e,T), ^^'^^ 

where r = is a slow variable. 

Let us substitute (13.21) into system (13.11) and group the coefficients of the 
same powers of the small parameter e. The equation for Y takes obviously 
the form 

dlY + Y = 0. 
The general solution of this equation is of the form 



Y{e, r) = k{T) exp(z^) + A;(r) exp{-te) 

with arbitrary function A;(r) to be chosen later on. 
The main term xq is determined from the equation 

Q 

degXo + fidexo + (A^ -Y)xo = f cos -, 
which can be rewritten in the form 



^ss^o + 4/i dgXQ + (g — 2r cos 2s)xq = 4/ cos (^s — -a j (3.3) 



with 

9 + a 

s = a = arg/c, 

q = 4A^, r = 2\k\, 

cf. (EH). 
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Let (fi, (p2 be a fundamental system of solutions of the homogeneous equa- 
tion corresponding to f l3.3p . Then any solution of f l3.3p can be represented in 
the form 

J so Ms')Ms') - Ms')Ms') V 2/ 

the denominator being the Wronsky determinant of the linearly independent 
system {v3i,(/?2}- 

The change of dependent variables xq = u exp(— 2/is) reduces the homo- 
geneous equation corresponding to (13. 3p 

dl^XQ + 4/i dsXQ + {q-2r cos 2s)xo = (3.4) 

to the Mathieu equation with slowly varying coefficient r. More precisely, we 
get 

n" + (g - 2i? cos 2s)n = 0, (3.5) 

where Q = q — ji^ and R = R{e^s). 

It is known that the general solution of the Mathieu equation changes 
within a period as ^1,2(5 -l- 2n) = exp(27rAi^2)wi,2(s), where Xi,2{(l,^) sue 
characteristic indices of the (homogeneous) Mathieu equation, see [10]. For 
q and r such that both the indices are purely imaginary, the solution of the 
Mathieu equation is bounded. If either of the characteristic indices is real 
and positive then the solution has exponential growth. 

Depending on the values q and r the solution can grow exponentially or 
remain bounded, see [12]. In Fig. Elthe dependence upon the characteristic 
index 3?Ai > on the parameters Q and R is shown. 

A section of the surface Xi{Q,R) for a fixed value Q and a line corre- 
sponding to dissipation are presented in Fig. HI where Q = 36, /i = 0.1. 

In the case under consideration the coefficient r = r(r) changes slowly. 
In each 2-7r -interval of the variable s the multiplicator of the solution of (13. 4p 
has the form 

/ci,2 = exp(27r(Ai,2(g, r(r)) - 2fi)). 

Consider the properties of solutions of (13. 4p for the fixed value g = 36, see 
Fig. m For 3f?A — 2/i < the multiplicators are less than 1. Hence it 
follows that the modulus of the solution decreases on a sequence of intervals. 
If 3fJA — 2/i > then the solution growths exponentially on a sequence of 
intervals. 

In Fig. |5]the integral index 



A{t)= f {m,{Q,R{r'))-2i2)dT' 
Jo 
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Fig. 3: The dependence of the real part of characteristic index 3?Ai of the 
Mathieu equation f l3.5p upon parameters Q and R. This figure is obtained by 
numerical simulations for different values of parameters, where Q G [25, 49] 
and R G [0, 64]. The grid step for numerical simulations equals 0.01 for both 
parameters. 
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Fig. 4: A section of the surface of the real part of characteristic index for 
the Mathieu equation for a fixed Q = 36. The dashed line is related to the 
dissipation in equation (13. 4|) . The coefficient of dissipation equals 0.1 for 
numerical simulations. 

is shown depending on time. The solution of equation (13. 4p does not depend 
in the main on the initial data as long as yl < 0. For equation (12. ip this 
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Fig. 5: The dependence of the integral index A upon 6. 

means that when studying the general solution one can restrict oneself to the 
particular solution of zero initial data to the inhomogeneous equation and 
neglect the solution of the homogeneous equation as long as the condition 
/I < is fulfilled. 

The correction yi is determined from the differential equation 

d'eeyi + yi = -2dlY + 6xl. 

We look for a solution yi of the form 



y^{e, t) = £(r) exp{te) + £{t) exp(-ze), 
satisfying the equation 

1 



doi = drk exp(— 2z6') — drk + — 5x1 exp(— z6'). 

If q and r are in the domain where the characteristic indices of the Mathieu 
equation are purely imaginary and 

drk = — lim — / xlexp(—i'd)d'd, (3.6) 
T^oo T Jo 

then 

yi = o{e) 

for 6* — 7- oo. The averaging over the fast variable on the right-hand side of 
(13.61) determines the derivative of k as function of slow time . The evolution 
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of Xq over the fast variable is determined by the inhomogeneous equation 
fl2.ip . The dependence of Xq upon the slow variable is explained by the 
slow perturbation of the coefficients of the Mathieu equation. The explicit 
dependence of xq upon the slow variable has not so far been defined. However, 
one can verify (13. 6p by means of numerical simulation of the initial system 

The numerical solution of system ( II. ip allows one to evaluate separately 
the left-hand side and the right-hand side of (13.61) and to determine the 
relative residual. We construct a numerical solution for the parameter values 
e = 0.2, f] = 1, u = 3, A = 1 and u = 0.1. 

The system (II. ip is solved by the Runge-Kutta method of 4 th order. As 
result we get a numerical solution Xnum, Unum- We then divide the entire 
integration interval into subintervals of length 27r and on each subinterval 
we compute the Fourier coefficients of sin 9 and cos 9 for the function eynum- 
From the data obtained in this way we evaluate the difference derivative k'^. 

The averaging operator on the right-hand side of (13. 6p is written in the 
form convenient for analytical computations. However, in numerical data 
the dependencies upon fast and slow variables are not separated from each 
other. Hence, it is not possible to directly apply formula (13. 6p . Instead the 
averaging operator over the fast variable of (13.61) is replaced by the averaging 
over the interval of length 0(£~^) with centre at the point Tj = e'^ti. In this 
way we get the values of the right-hand side at the points Xj. 

In Fig. El one compares the derivatives in r on the left-hand side of 
equation (13. 6p . which are evaluated numerically, and the integral on the right- 
hand side. 

We thus conclude that the substitution of the numerical solution for the 
genuine solution of (13. 6p leads to an inessential residual. 



4 Analysis of asymptotic solution for A ^ 1 

In this section we carry out the analysis of the behaviour of the main terms 
of representation (13. 2p under the assumption that X = u/ f2 ^ 1. This 
assumption corresponds to the case of forced combinative scattering, where 
A attains the value 10^, see [3], [1]. 

Equation (12.11) for the main term xo can be obviously rewritten in the 
form 

o,,Xq + — dsXo + (^4 - — cos 2s)xo = ^ 2)' 
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Fig. 6: The relative error jS"; — ^rl/IS"/! of equation (13.61) . Here Si and Sr are 
numerical evaluations of the left-hand and right-hand sides of (13. 6p . respec- 
tively. 



The asymptotics of the particular solution for large values A is 

2/ cos(s + a/2) 



Xo 



2 — A~2r cos 2s 



x-' + o{x- 



(4.1) 



Our next objective is to treat equation (13.61) for A ^ 1. To this end we 
substitute (14. ip into equation (13.60 and evaluate the integral explicitly. As 
result we get 



dr' 



OH 

dK'2 
dH 



(4.2) 



where r' = A ^t. This is a Hamiltonian system with Hamiltonian 



H 



\K\ 



2tt[2^ Ki(1 - |K|2) - {\K\^ - Ki)^l - \K\^' 



The system (14. 2 p is easily seen to have a stable equilibrium point of type 
'centre' 

K2 = 0, 
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Fig. 7: A phase portrait of solutions of (14. 2 p in a neighborhood of 'centre' 
Ki = 1/V2, K2 = 0. This fi gure is obtained numerically with diverse initial 
data for system (14. 2p . 



see Fig. [71 

The amplitude of oscillations can be found from the system 



d 



Ho, 
0. 



(4.3) 



In Fig. [H] the dependence of the largest r+ (dashed curve) and smallest r_ 
(solid curve) values of solutions of (14. 3 p upon the Hamiltonian is shown. The 
vertical line corresponds to Hq = l/(47r). This allows one to evaluate the 
maximal value of the amplitude of oscillation envelope through 



max{|r_|, |r_|_|}, 



where 



The period is determined by 



(4.4) 



dKi 

Wo 



(4.5) 
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Fig. 8: Dependence of r+ (dashed) and r_ (solid) evaluated by fl4.4p on Hq. 
The vertical line corresponds to the value Hq = l/(4-7r). An intersection of 
these lines allows one to estimate the maximal value of envelope amplitude. 

where the integration is over the cycle c in the complex plane of the variable 
K = Ki + given by 
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In Fig. [9] one sees the results of numerical simulation for the solution of 
system (11. ip with initial data corresponding to the stationary solution 

^(0) = ^, 1/(0) 



A2' ''^^ eV2' 
x'(0) = 0, y'{0) = 0. 

The analytical result obtained in the paper agrees well to the numerical 
simulation. For the values of parameters Ki and K2, which are used in the 
numerical simulation, the Hamiltonian value proves to be Hq = l/(47r). In 
this case the amplitude of pulsations is evaluated by X'^e~^. It is shown in 
Fig. [m The horizontal line corresponds to A^e"^. 

The explicit value for the period of envelope oscillations can be evaluated 
from (14. 5p . Our numerical simulations with parameters A = l,u = 3,fl = 
1, /i = 0.1, e = 0.2 give the value T = 5264.76 in the variable r. 

In Fig. [TU] the pulsations of the y -component are shown as well as the 
envelope superposed on them. The envelope is computed numerically as 
solution of system (14. 2 p with e = 0.2, \ = 3, A = 1. All these numerical 
evaluations are realized by the Adams method with the precision 10"^". 
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Fig. 9: The constant envelope function and oscillations of the component y. 
The component y is obtained by numerical simulations for the solution of sys- 
tem (11. ip . The constant envelope function (dashed curve) equals ±A^e~^/A/2 
which is related to the stationary solution Ki = 1/ \/2, = of Hamiltonian 
system (H^ . 
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Fig. 10: The numerically evaluated y -component of the solution and the 
envelope. The envelope function is evaluated as the solution of fl4.2l) . 

System (14. 2 p gives an approximation for determining the behaviour of 
functions Ki and K2- In Fig. [TT]the solid curve corresponds to the envelope 
of the y -component of the solution. This curve is found from a numerical so- 
lution of system (II. ip . The dashed curve \J Kf + gives an approximation 
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Fig. 11: Numerically evaluated envelope function of ?/ -component of the 
solution (solid curve) and the approximation of the solution envelope as so- 
lution of system (14. 2 p (dashed curve). The vertical line corresponds to the 
value of the period T = 5264.76. It is determined numerically by (14. 5p . The 
horizontal line is an estimate of maximal value A^e"^ of envelope amplitude. 

of the solution envelope as solution of system (14. 2p . A point Ki = 0, K2 = 
is an unstable node equilibrium point. Initial data for system (14. 2 p corre- 
spond to the linear resonance of ?/ -component of the solution for (II. ip in 
the initial stage, see Fig. [T] In Fig. [TT]the vertical line shows the value of 
the period T = 5264.76. It is determined numerically by expression (14.51) . 
The horizontal line is an estimate of the maximal value A^e:"^ of envelope 
amplitude. 

5 Conclusion 

It is shown in the paper that for the large time t ~ with e being a 
coupling parameter the behaviour in the main of a resonantly perturbed 
system of coupled oscillators is described by a Hamiltonian system. This 
system is studied in the physically interesting case when the quotient of 
oscillator frequencies w/i? is large enough. Certain pulsations prove to occur 
in the solution of the initial system of coupled oscillators. The Hamiltonian 
system obtained allows one to evaluate the period and amplitude of these 
pulsations. 
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